unit CannyEdge;


interface

uses
 Classes, SysUtils, Math, GR32, GR32_Filters, GR32_VectorMaps, GR32_OrdinalMaps,
 GR32_Resamplers, GR32_Containers, GR32_FloatMap, GR32_LowLevel;

const
  MAX_KERNEL_SIZE = 12;

type
  TKernelSize = 1..MAX_KERNEL_SIZE;

  TKernel = record
    Size: TKernelSize;
    Weights: array[-MAX_KERNEL_SIZE..MAX_KERNEL_SIZE] of Single;  
  end;
type
 TGrayMethod = (gmIntensity, gmUniform, gmRed, gmGreen, gmBlue);
 TCannyEdge = Class
 Private
  fKarnel :Integer;
  fLowThresh, fHighThresh: Single;
  fScale, fOffset : Integer;
  fAutoThresh     : Boolean;
  procedure MakeGaussianKernel(var K: TKernel; Size: Integer);
 Public
  procedure CannyEdgeDetect(Src, Dst: TBitmap32;
    GrayMethod: TGrayMethod = gmIntensity);
  procedure CannyEdgeDetectEx(Src, Dst: TBitmap32);
  procedure GaussianFilter(Src: TIntegerMap; Size: integer);

  Constructor Create;

  published
    property KanelSize: integer read fKarnel write fKarnel;
    property LowTresh : Single  read fLowThresh write fLowThresh;
    property HighTresh: Single  read fHighThresh write  fHighThresh;
    property Scale    : Integer read fScale write fScale;
    property Offset   : Integer read fOffset write  fOffset;
    property AutoThresh: Boolean read fAutoThresh write fAutoThresh;
end;



implementation

uses
   Contnrs;

Constructor TCannyEdge.Create;
begin
  Inherited Create;
  fKarnel     :=    2;
  fLowThresh  :=    8;
  fHighThresh :=    16;
  fScale      :=    16;
  fOffset     :=    128;
  fAutoThresh := False;

end;
function MeanValue(Map: TIntegerMap): Single;  
var  
  I, Count: Integer;  
  Mean: Single;  
begin  
  Mean := 0.0;  
  Count := 0;  
  for I := 0 to Map.Width * Map.Height - 1 do  
  begin  
    if Map.Bits[I] > 0 then  
    begin  
      Mean := Mean + Map.Bits[I];  
      Inc(Count);  
    end;  
  end;  
  Mean := Mean / Count;  
  Result := 0;  
  Count := 0;  
  for i := 0 to Map.Width * Map.Height - 1 do  
  begin  
    if Map.Bits[I] > Mean then  
    begin  
      Result := Result + Map.Bits[I];  
      Inc(Count);  
    end;  
  end;  
  Result := Result / Count;  
end;  
function GrayValue(Color32: TColor32; GrayMethod: TGrayMethod): Integer;  
begin  
  case GrayMethod of  
    gmUniform:  
      Result := ((Color32 and $FF0000) shr 16 + (Color32 and $FF00) shr 8 + (Color32 and $FF)) div 3;  
    gmRed:  
      Result := (Color32 and $FF0000) shr 16;  
    gmGreen:  
      Result := (Color32 and $FF00) shr 8;  
    gmBlue:  
      Result := Color32 and $FF;  
  else  
    Result := Intensity(Color32);  
  end;  
end;  
function Sector(Theta: Double): Integer;  
const  
  C = 180 / Pi;  
begin  
  Theta := Theta * C;  
  Theta := Round(Theta + 270) mod 360;  
  Result := 0;  
  if ((theta >= 337.5) or (theta < 22.5) or ((theta >= 157.5) and (theta < 202.5))) then  
    Result := 0  
  else if (((theta >= 22.5) and (theta < 67.5)) or ((theta >= 202.5) and (theta < 247.5))) then  
    Result := 1  
  else if (((theta >= 67.5) and (theta < 112.5)) or ((theta >= 247.5) and (theta < 292.5))) then  
    Result := 2  
  else if (((theta >= 112.5) and (theta < 157.5)) or ((theta >= 292.5) and (theta < 337.5))) then  
    Result := 3  
end;  
function Suppress(Map: TIntegerMap; Sector, I, J: Integer; LowThresh: Single): Integer;  
var  
  Tmp: Integer;  
begin  
  Tmp := Map[I, J];  
  Result := 0;  
  if Tmp < LowThresh then Exit;  
  if (Sector = 0) then  
  begin  
    if not ((Map[I + 1, J] > Tmp) or (Map[I - 1, J] >= Tmp)) then  
      Result := Tmp;  
  end  
  else if (Sector = 1) then  
  begin  
    if not ((Map[I + 1, J + 1] > Tmp) or (Map[I - 1, J - 1] >= Tmp)) then  
      Result := Tmp  
  end  
  else if (Sector = 2) then  
  begin  
    if not ((Map[I, J + 1] > Tmp) or (Map[I, J - 1] >= Tmp)) then  
      Result := tmp;  
  end  
  else if (Sector = 3) then  
  begin  
    if not ((Map[I + 1, J - 1] > Tmp) or (Map[I - 1, J + 1] >= Tmp)) then  
      Result := Tmp;  
  end;  
end;
procedure Track(Src: TIntegerMap; Marked: TBooleanMap; Stack: TStack;  
  Thresh: Single; I, J: Integer);  
var  
  A, B: Integer;  
begin

  Stack.Push(Pointer(I));
  Stack.Push(Pointer(J));
  while (Stack.Count > 0) do  
  begin  

    B := Integer(Stack.Pop);
    A := Integer(Stack.Pop);  

    if (Marked[A, B] = False) then
    begin  

      if (A > 0) and (B > 0) then
        if (Src[A - 1, B - 1] > Thresh) then
        begin  
          Stack.Push(Pointer(A - 1));  
          Stack.Push(Pointer(B - 1));  
        end;  
      if (B > 0) then  
        if (Src[A, B - 1] > Thresh) then  
        begin  
          Stack.Push(Pointer(A));  
          Stack.Push(Pointer(B - 1));  
        end;  
      if (A < Src.Width - 1) and (B > 0) then  
        if (Src[A + 1, B - 1] > Thresh) then  
        begin  
          Stack.Push(Pointer(A + 1));  
          Stack.Push(Pointer(B - 1));  
        end;  
      if (A > 0) then  
        if (Src[A - 1, B] > Thresh) then  
        begin  
          Stack.Push(Pointer(A - 1));  
          Stack.Push(Pointer(B));  
        end;  
      if (A < Src.Width - 1) then  
        if (Src[A + 1, B] > Thresh) then  
        begin  
          Stack.Push(Pointer(A + 1));  
          Stack.Push(Pointer(B));  
        end;  
      if (A > 0) and (B < Src.Height - 1) then  
        if (Src[A - 1, B + 1] > Thresh) then  
        begin  
          Stack.Push(Pointer(A - 1));  
          Stack.Push(Pointer(B + 1));  
        end;  
      if (B < Src.Height - 1) then  
        if (Src[A, B + 1] > Thresh) then  
        begin  
          Stack.Push(Pointer(A));  
          Stack.Push(Pointer(B + 1));  
        end;  
      if (A < Src.Width - 1) and (B < Src.Height - 1) then  
        if (Src[A + 1, B + 1] > Thresh) then  
        begin  
          Stack.Push(Pointer(A + 1));  
          Stack.Push(Pointer(B + 1));  
        end;  
      Marked[A, B] := True;  
    end;  
  end;  
end;  

procedure ApplyTrack(Src, Dst: TIntegerMap; LowThresh, HighThresh: Single);  
var  
  Marked: TBooleanMap;  
  Stack: TStack;  
  I, J: Integer;  
begin  
  Marked := TBooleanMap.Create;  
  Stack := TStack.Create;  
  try  
    Marked.SetSize(Src.Width, Src.Height);  
    Marked.Clear(0);  

    for I := 0 to Src.Width - 2 do
      for J := 0 to Src.Height - 2 do

        if ((Src[I, J] > HighThresh) and (Marked[I, J] = False)) then  
          Track(Src, Marked, Stack, LowThresh, I, J);

    for I := 0 to Src.Width - 1 do  
      for J := 0 to Src.Height - 1 do  
      begin  
        if (Marked[I, J] = False) then  
          Dst[I, J] := 0  
        else  
          Dst[I, J] := Src[I, J];  
      end;  
  finally  
    if Assigned(Stack) then Stack.Free;  
    if Assigned(Marked) then Marked.Free;  
  end;  
end;
Procedure TCannyEdge.CannyEdgeDetect(Src, Dst: TBitmap32;
    GrayMethod: TGrayMethod = gmIntensity);
var
  Temp: TIntegerMap;
  P, Q, M: TIntegerMap;
  ThetaMap: TFloatMap;
  SectorMap, NMSMap: TIntegerMap;
  I, J, R: Integer;
  W, H: Integer;
  Median: Double;
begin
  W := Src.Width;
  H := Src.Height;
  if W * H = 0 then Exit;
  //Setup local variables;
  Temp := TIntegerMap.Create;
  P := TIntegerMap.Create;
  Q := TIntegerMap.Create;
  M := TIntegerMap.Create;
  ThetaMap := TFloatMap.Create;
  try
    Temp.SetSizeFrom(Src);
    P.SetSizeFrom(Src);
    Q.SetSizeFrom(Src);  
    M.SetSizeFrom(Src);  
    ThetaMap.SetSizeFrom(Src);  

    for I := 0 to Src.Width - 1 do  
      for J := 0 to Src.Height - 1 do  
        Temp[I, J] := GrayValue(Src[I, J], GrayMethod);  

    GaussianFilter(Temp, fKarnel);

    for I := 1 to W - 2 do  
      for J := 1 to H - 2 do  
      begin  
        P[I, J] := (Temp[I, J + 1] - Temp[I, J] +  
          Temp[I + 1, J + 1] - Temp[I + 1, J]) div 2;  
        Q[I, J] := (Temp[I, J] - Temp[I + 1, J] +  
          Temp[I, J + 1] - Temp[I + 1, J + 1]) div 2;  
        M[I, J] := Round(Sqrt(P[I, J] * P[I, J] + Q[I, J] * Q[I, J]));  
        ThetaMap[I, J] := ArcTan2(Q[I, J], P[I, J]);  
      end;  

    if fAutoThresh then
    begin  
      Median := MeanValue(M);  
      fLowThresh := Round(0.5 * Median);
      fHighThresh := Round(Median);
    end;  

    SectorMap := P;  
    SectorMap.Clear(0);  
    for I := 0 to W - 2 do  
      for J := 0 to H - 2 do  
        SectorMap[I, J] := Sector(ThetaMap[I, J]);

    NMSMap := Q;
    NMSMap.Clear(0);  
    for I := 1 to W - 2 do  
      for J := 1 to H - 2 do  
        NMSMap[I, J] := Suppress(M, SectorMap[I, J], I, J, fLowThresh);

    ApplyTrack(NMSMap, M, fLowThresh, fHighThresh);
    Dst.SetSizeFrom(Src);

    for I := 0 to W - 1 do  
      for J := 0 to H - 1 do  
      begin  
        R := M[I, J];
        R := R * fScale;
        if R <> 0 then R := R + fOffset;
        if (R > 255) then R := 255;  
        if (R < 0) then R := 0;  
        Dst[I, J] :=  (R shl 16) or (R shl 8) or R;
      end;  
    Dst.Changed;  
  finally  
    if Assigned(Temp) then Temp.Free;  
    if Assigned(P) then P.Free;  
    if Assigned(Q) then Q.Free;  
    if Assigned(M) then M.Free;  
    if Assigned(ThetaMap) then ThetaMap.Free;  
  end;
end;

Procedure TCannyEdge.CannyEdgeDetectEx(Src, Dst: TBitmap32);
var  
  Temp1, Temp2, Temp3: TBitmap32;
  I, J: Integer;  
  C1, C2, C3: LongInt;  
begin  
  Temp1 := TBitmap32.Create;
  Temp2 := TBitmap32.Create;
  Temp3 := TBitmap32.Create;
  try  
    Temp1.SetSizeFrom(Src);  
    CannyEdgeDetect(Src, Temp1, gmRed);
    Temp2.SetSizeFrom(Src);
    CannyEdgeDetect(Src, Temp2, gmGreen);
    Temp3.SetSizeFrom(Src);  
    CannyEdgeDetect(Src, Temp3, gmBlue);
    Dst.SetSizeFrom(Src);  
    // Combine the result  
    for I := 0 to Src.Width - 1 do  
      for J := 0 to Src.Height - 1 do  
      begin  
        C1 := Temp1[I, J] ;  
        C2 := Temp2[I, J] ;  
        C3 := Temp3[I, J] ;  
        Dst[I, J] := C1 or C2 or C3;  
      end;  
    Dst.Changed;  
  finally  
    Temp1.Free;  
    Temp2.Free;  
    Temp3.Free;  
  end;  
end;
Procedure TCannyEdge.MakeGaussianKernel(var K: TKernel; Size: Integer);
var  
  I: Integer;  
  Temp: Double;  
begin  
  K.Size := Size;  

  K.Weights[0] := 1.0;

  for I := 1 to Size do  
  begin  
    Temp := 3.0 * I / Size;  
    Temp := Exp(-Temp * Temp / 2.0);  
    K.Weights[I] := Temp;  
    K.Weights[-I] := Temp;  
  end;  

  Temp := 0;
  for I := -K.Size to K.Size do  
  begin  
    Temp := Temp + K.Weights[I];  
  end;  
  for I := -K.Size to K.Size do  
  begin  
    K.Weights[I] := K.Weights[I] / Temp;  
  end;  
end;

Procedure TCannyEdge.GaussianFilter(Src: TIntegerMap; Size: Integer);
var  
  Col, Row: Integer;  
  K: TKernel;  
  function ApplyGaussianH(Col, Row: Integer): Single;  
  var  
    I: Integer;  
  begin  
    Result := 0;  
    for I := -K.Size to K.Size do  
    begin  
      Result := Result + Src[Clamp(Col + I, 0, Src.Width - 1), Row] *  
        K.Weights[I]  
    end;  
  end;  
  function ApplyGaussianV(Col, Row: Integer): Single;  
  var  
    I: Integer;  
  begin  
    Result := 0;  
    for I := -K.Size to K.Size do  
    begin  
      Result := Result + Src[Col, Clamp(Row + I, 0, Src.Height - 1)] *  
        K.Weights[I]  
    end;  
  end;  
begin  
  MakeGaussianKernel(K, Size);  
  for Row := 0 to Src.Height - 1 do  
    for Col := 0 to Src.Width - 1 do  
      Src[Col, Row] := Round(ApplyGaussianH(Col, Row));  
  for Row := 0 to Src.Height - 1 do  
    for Col := 0 to Src.Width - 1 do  
      Src[Col, Row] := Round(ApplyGaussianV(Col, Row));  
end;
end.
